function nicfps_unwarp2,inarr,area_weight=area_weight,weights=weights, $
    _extra=extra, nan_border=nanbord

;+
; NAME:
;     NICFPS_UNWARP (Revision 2)
; PURPOSE:
;     Remove geometric distortion from one or more NIC-FPS images.
; EXPLANATION:
;     Uses the IDL function POLY_2D and a 2D polynomial fit to the NIC-FPS
;     distortion map to generate new, un-distorted image(s) from one or more
;     NIC-FPS images.  By default,
;     this method does NOT maintain flux; normally, that is achieved by flat-
;     fielding.  This routine simply re-samples (and interpolates) the values
;     in the input image(s) at the appropriate locations to determine the
;     values in each pixel of the output image(s).  An optional keyword can
;     be used to change this behavior in order to maintain flux.
; CALLING SEQUENCE
;     outarr = NICFPS_UNWARP ( inarr [, /nan_border, /area_weight, $
;                            weights=weights ] )
; INPUT:
;     inarr   -- A 2D image or a 3D array containing a stack of 2D NIC-FPS
;                images.  These should have the same geometry as raw NIC-FPS
;                images, i.e. 1024x1024 images with North up, East left.
;                (Reference pixels should not be cropped.)
; OUTPUT:
;     outarr  -- An array of the same dimensions and type as INARR, containing
;                the input images warped to remove geometric distortion.  If
;                the keyword NAN_BORDER is set, OUTARR is of type FLOAT.
; OPTIONAL KEYWORDS:
;     nan_border -- Setting this keyword has two effects: 1) reference pixels
;                   are set to NaN (!values.f_nan), and 2) regions in each
;                   new image which correspond to locations outside INARR
;                   are set to NaN.  (Setting these borders to NaN allows
;                   this border region to be ignored in future averaging,
;                   median combining, etc.)
;     area_weight-- setting this keyword will scale the values of OUTARR such
;                   that all the flux in INARR is conserved.  Effectively
;                   multiplies each pixel value in OUTARR by the area (in
;                   pixels) in INARR to which that pixel corresponds.
;                   Usually, this correction is achieved by flat-fielding, in
;                   which case this keyword should NOT be set.
;     weights    -- If the AREA_WEIGHT keyword is set, this keyword can be
;                   set to the name of a variable which on output will contain
;                   the weights applied to OUTARR.  This will be a FLOAT array
;                   of the same dimensions as OUTARR (although each 2D slice
;             of WEIGHTS will be identical).
; EXAMPLE:
;     The following procedure will remove distortion from three NIC-FPS images
;     which have been read into the variables im1,im2,im3:
;
;     IDL> imstack = [ [[im1]], [[im2]], [[im3]] ]  ;stack images in 3D array
;     IDL> unwarped = nicfps_unwarp(imstack, /nan_border)
;
;     the variable unwarped now contains the un-distorted image, with edges
;     (reference pixels and regions which map outside of input images)
;     set to NaN (Not a Number).
; PROCEDURE:
;     Locations (x',y') in each new image correspond to locations (x,y) in
;     the corresponding input image according to the 2D polynomial equations:
;     x'=sum(coeff_x[i,j]*x^j*y^i), y'=sum(coeff_y[i,j]*x^j*y^i),
;     with coefficients defined below.  The value at (x',y') is then
;     set by POLY_2D to the value interpolated for (x,y) in the input image.
;
;     /NaN_BORDER sets MISSING=!values.f_nan in the call to POLY_2D.
;     /AREA_WEIGHT computes the 1st partial derivatives of the 2D polynomial
;     functions; these give the delta(x) and delta(y) (in input pixels) of
;     each output pixel.  delta(x)*delta(y) is taken to be the area of the
;     region to which each new pixel corresponds, and is the factor used to
;     adjust the value of each new pixel.
; MODIFICATION HISTORY:
;     Written by Nathaniel Cunningham, Univ. of Colorado, June 2005, using
;         polynomial coefficients determined by A. Bondarenko.
;-

;polynomial coefficients for poly_2d.
;x'=sum(coeff_x[i,j]*x^j*y^i), y'=sum(coeff_y[i,j]*x^j*y^i)

   ;NEW COEFFICIENTS
   coeff_x = [[-7.6423547, 0.018567389, -1.6360355e-005, 1.6949539e-009], $
              [1.0279830, -2.9072323e-005, 2.6017212e-008, 0.0000000000], $
              [-4.4690867e-005, -1.0943165e-009, 0.0000000, 0.000000000], $
              [2.8775651e-008, 0.00000, 0.00000, 0.00000]]


   coeff_y = [[-7.4005260,  1.0296221, -3.8010270e-005, 2.1732379e-008], $
              [0.014171051, -3.0130551e-005, 7.4258827e-010, 0.0000000], $
              [-1.3542434e-005, 2.7194525e-008, 0.00000000, 0.00000000], $
              [-1.2327870e-009, 0.000000000000, 0.000000000, 0.0000000]]



  ;OLD COEFFICIENTS
  ; coeff_x = [[-7.1042072, 0.017056966, -1.3799720e-005, -5.5723316e-011], $
  ;             [1.0288530, -3.1588962e-005, 2.5825762e-008, 0.0000000000], $
  ;             [-4.3730946e-005, -1.3832661e-009, 0.00000000, 0.00000000], $
  ;             [2.8547903e-008, 0.0000000000, 0.0000000000, 0.0000000000]]

  ;  coeff_y = [[-10.542904, 1.0347672, -4.5872982e-005, 2.4582113e-008], $
  ;             [0.017423202, -2.6262299e-005, -2.2224779e-010, 0.00000], $
  ;             [-1.5369865e-005, 2.4949451e-008, 0.0000000, 0.00000000], $
  ;             [-2.5236412e-010, 0.0000000000, 0.0000000000, 0.0000000]]

;check input array dimensions and type
ndims=size(inarr,/n_dimensions)
if ndims lt 2 or ndims gt 3 then begin
    print, 'INARR must be a 2-D image or a 3-D stack of images'
    return,-1
endif
help,inarr
dims=size(inarr,/dimensions)
if ndims eq 3 then nims=dims[2] else nims=1     ;how many images?

arrtype=size(inarr,/type)     ;what datatype?

;create output array of matching dimensions and type
outarr=make_array(dims,type=arrtype,/nozero)

;un-warp each image in the stack
if keyword_set(nanbord) then begin
    newarr=make_array(dims,value=!values.f_nan)
    newarr[4:1019,4:1019,*]=inarr[4:1019,4:1019,*]
    for i=0,nims-1 do outarr[0,0,i]=poly_2d(newarr[*,*,i], $
       coeff_x,coeff_y,missing=!values.f_nan,cubic=-0.5,_extra=extra)
endif else begin
    for i=0,nims-1 do outarr[0,0,i]=poly_2d(inarr[*,*,i], $
       coeff_x,coeff_y,cubic=-0.5,_extra=extra)
endelse

if keyword_set(area_weight) then begin
    ncols=dims[0]
    nrows=dims[1]
    cdims=size(coeff_x,/dimensions)
    nccol=cdims[0]
    ncrow=cdims[1]
    xpos=rebin(findgen(ncols),ncols,nrows)
    ypos=rebin(reform(findgen(nrows),1,nrows),ncols,nrows,nccol)

    ywid=(xwid=fltarr(ncols,nrows))
    for icol=0,3 do begin
       for irow=1,3 do begin
         xwid += coeff_x[icol,irow]*float(irow)*xpos^(irow-1)*ypos^icol
         ywid += coeff_y[irow,icol]*xpos^icol*float(irow)*ypos^(irow-1)
       endfor
    endfor
    weights=rebin(xwid*ywid,ncols,nrows,nims)
    outarr=outarr*weights
endif

return,outarr

end
